# -*- coding: utf-8 -*-
"""
@author: asus

用于解决包含 31 个城市的 TSP 问题。

问题的规划模型表述如下：
"""

# In[0]:
# 导入数值计算工具包
import numpy as np

# 导入 Pyomo 环境和求解器调用工具
from pyomo.environ import *
from pyomo.opt import SolverFactory

# 导入自定义的数据处理过程
from data_process import load_tsp_data, get_adj_mat

# In[1]:
# 读取数据并生成边权矩阵
adj_mat = get_adj_mat(
    load_tsp_data())

# N 为边权矩阵的长度
N = len(adj_mat)

# In[2]:
# 创建规划模型对象
model = ConcreteModel()

# In[3]:
# 创建规划变量
model.i = RangeSet(1, N-1)
model.j = RangeSet(2, N)

model.ui = RangeSet(1, N)

model.x = Var(model.i, model.j, domain=Binary)
model.u = Var(model.ui, domain=NonNegativeIntegers)

# In[4]:
# 创建规划目标

model.obj = Objective(
    expr=(
        sum(
            sum(
                model.x[i, j] * adj_mat[i-1][j-2] for i in model.i
                ) for j in model.j)),
    sense=minimize)

# In[5]:
# 创建模型约束
def con_x_rule(model, i):
    return sum(model.x[i, j] for j in model.j) == 1

def con_y_rule(model, j):
    return sum(model.x[i, j] for i in model.i) == 1

model.con_x = Constraint(model.i, rule=con_x_rule)
model.con_y = Constraint(model.j, rule=con_y_rule)

def con_u_rule(model, i, j):
    if (i != j):
        return (model.u[i] - model.u[j] + N * model.x[i, j]) <= N-1
    elif (i == j):
        return model.x[i, j] == 0

model.con_u = Constraint(
    model.i, model.j, rule=con_u_rule)

# In[6]:
# 求解模型

opt = SolverFactory('scip') # 调用 SCIP
solver_options = opt.options

results = opt.solve(model) # 求解模型

print(results)

# In[7]:

model.pprint()

# In[7]:
# 打印输出结果

for i in range(1, N):
    for j in range(2, N+1):
        if value(model.x[i, j]) == 1:
            print(i, '-->', j)

# In[8]:
for i in range(1, 32):
    print(value(model.u[i]))